Spectra Story

Author

Adam Kemberling

Published

October 18, 2024

Code
# Libraries
library(tidyverse)
library(gmRi)
library(rnaturalearth)
library(scales)
library(sf)
library(gt)
library(patchwork)
library(lmerTest)
library(emmeans)
library(merTools)
library(tidyquant)
library(ggeffects)
library(performance)
library(gtsummary)
library(gt)
library(sizeSpectra)
library(ggdist)
library(pander)
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflict_prefer("filter", "dplyr")
conflicted::conflicts_prefer(lmerTest::lmer)

# Processing functions
source(here::here("Code/R/Processing_Functions.R"))

# ggplot theme
theme_set(
  theme_gmri(
    rect = element_rect(fill = "white", color = NA), 
    strip.text.y = element_text(angle  = 0),
    axis.text.x = element_text(size = 8),
    legend.position = "bottom"))



# vectors for factor levels
area_levels <- c("GoM", "GB", "SNE", "MAB")
area_levels_long <- c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")
area_levels_all <- c("Northeast Shelf", area_levels)
area_levels_long_all <- c("Northeast Shelf", area_levels_long)

# table to join for swapping shorthand for long-hand names
area_df <- data.frame(
  area = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "All"),
  survey_area = c("SS", "GoM", "GB", "SNE", "MAB", "Northeast Shelf"),
  area_titles = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"))


# Degree symbol
deg_c <- "\u00b0C"



# Function to get the Predictions
# Drop effect fits that are non-significant  ###
get_preds_and_trendsignif <- function(mod_x){
  modx_preds <- as.data.frame(
    # Model predictions
    ggpredict(
      mod_x, 
      terms = ~ est_year * survey_area * season) ) %>% 
    mutate(
      survey_area = factor(group, levels = area_levels),
      season = factor(facet, levels = c("Spring", "Fall")))
  
    # Just survey area and year
    modx_emtrend <- emtrends(
      object = mod_x,
      specs =  ~ survey_area*season,
      var = "est_year") %>% 
      as_tibble() %>% 
      mutate(
        zero = 0,
        non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))
  
    # Preds with signif
    modx_preds %>% left_join(select(modx_emtrend, survey_area, season, non_zero))
  
}

Storycrafting - Northeast Shelf Spectra

I think we need to step back a bit and find the main story again–which is about how size structure of the NES fish community has changed over time, whether that varies by region, and what covariates are associated with those changes.
The Aug draft you developed really focused on comparing regression models to consider drivers based on the fish community size spectrum as viewed from the perspective of the Wigley species set vs. the full species set.
However, I didn’t see figures of those spectra in that draft (maybe there are somewhere else and I’ve lost track). The recent slide deck really focuses on seasonal differences in the size spectra and influence of temperature. - KMills

This document will be split into two sections (the second potentially moving to its own doc).

Part 2: Relationships to External Factors

The aim of this section is to provide context and potentially attribution/correlations that may help explain trends/changes observed in the above section.

Code
# # Model Dataset for Wigley Species Subset
wigley_bmspectra_model_df <- read_csv(here::here("Data/model_ready/wigley_community_bmspectra_mod.csv"))

# Use one for plots
wigley_bmspectra_model_df <- wigley_bmspectra_model_df %>% 
  mutate(survey_area = case_when(
    survey_area == "GoM" ~ "Gulf of Maine",
    survey_area == "GB" ~ "Georges Bank",
    survey_area == "SNE" ~ "Southern New England",
    survey_area == "MAB" ~ "Mid-Atlantic Bight"),
  survey_area = factor(survey_area, levels = area_levels_long),
  season = fct_rev(season))

Bottom Temperature Changes

Bottom temperatures are from the DuPontavice and Saba combined bottom temperature product. This data product has been clipped to the study areas and seasonally averaged to match the survey sampling seasons.

Code
# Model linear trends in time
btemp_mods <- lmer(
  formula = bot_temp ~ est_year * survey_area * season + (1|est_year),
  data = wigley_bmspectra_model_df)



# change the stupid labels around because
btemp_preds <- as.data.frame(
    # Model predictions
    ggpredict(
        btemp_mods, 
        terms = ~ est_year * survey_area * season) ) %>% 
    mutate(
        survey_area = factor(group, levels = area_levels_long),
        season = factor(facet, levels = c("Spring", "Fall")))
btemp_signif <- emtrends(
      object = btemp_mods,
      specs =  ~ survey_area*season,
      var = "est_year") %>% 
      as_tibble() %>% 
      mutate(
        zero = 0,
        non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))
btemp_trend_predictions <- left_join(btemp_preds, btemp_signif)  %>% 
  mutate(survey_area = group,
         survey_area = factor(survey_area, levels = area_levels_long),
         season = factor(season, levels = c("Spring", "Fall")))




# Plot temperatures
temp1 <- wigley_bmspectra_model_df %>% 
  ggplot() +
  geom_ribbon(
    data = filter(btemp_trend_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    aes(est_year, bot_temp, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, bot_temp, color = season, linetype = "5-Year Smooth"),
    size = 1, n = 5, ma_fun = SMA) +
  geom_line(
    data = filter(btemp_trend_predictions, non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_color_gmri() + 
  scale_fill_gmri() + 
  #scale_linetype_manual(values = c(1,1)) +
  facet_grid(survey_area~., scales = "free") +
  theme(legend.position = "right") +
  scale_y_continuous(labels = label_number(suffix = "\u00b0C")) +
  labs(y = "Bottom Temperature", x = "Year", color = "Season",
       fill = "Season", linetype = "Trend")



temp1 

Bottom temperatures have increased in all four regions (annual averages), but seasonally spring temperatures in SNE and MAB are stable.

Bottom Temperature Distributions

The distribution of seasonal bottom temperatures are such that Gulf of Maine experiences a much smaller difference in seasonal temperature swings than the other regions. In the Gulf of Maine the difference in bottom temperature from Spring to Fall can be smaller than the differences in spring bottom temperature along the whole Northeast Shelf.

Code
# Plot the distribution as a raincloud plot
wigley_bmspectra_model_df %>% 
  ggplot(aes(fct_rev(survey_area), bot_temp, color = season, fill = season)) +
   ggdist::stat_halfeye(
     alpha = 0.4,
     adjust = .5, 
     width = .6, 
    .width = 0, 
    justification = -.2, 
    point_colour = NA
  ) + 
  geom_boxplot(
    fill = "transparent",
    width = .15, 
    outlier.shape = NA
  ) +
  ## add dot plots from {ggdist} package
  ggdist::stat_dots(
    ## orientation to the left
    side = "left", size = 1,
    ## move geom to the left
    justification = 1.12, 
    ## adjust grouping (binning) of observations 
    binwidth = .1) + 
  scale_color_gmri() + 
  scale_fill_gmri() + 
  scale_y_continuous(labels = label_number(suffix = "\u00b0C")) +
  coord_flip() +
  theme(legend.position = "none") +
  labs(y = "Bottom Temperature", x = "", color = "Season")

The overlap in bottom temperature distributions is very similar to the overlap in mass spectra exponents. Gulf of Maine spectra and bottom temperatures share a good region of overlap, and the regions that experience larger seasonal fluctuations in BT see larger seasonal fluctuations in their size distributions.

Code
wigley_b_dist_plot <- wigley_bmspectra_model_df %>% 
  ggplot(aes(fct_rev(survey_area), b, color = season, fill = season)) +
  ggdist::stat_halfeye(
     alpha = 0.4,
     adjust = .5, 
     width = .6, 
    .width = 0, 
    justification = -.2, 
    point_colour = NA) + 
  geom_boxplot(
    fill = "transparent",
    width = .15, 
    outlier.shape = NA) +
  ## add dot plots from {ggdist} package
  ggdist::stat_dots(
    ## orientation to the left
    side = "left", size = 1,
    ## move geom to the left
    justification = 1.12, 
    ## adjust grouping (binning) of observations 
    binwidth = .005) + 
  scale_color_gmri() + 
  scale_fill_gmri() + 
  coord_flip() +
  labs(
    y = "ISD Exponent", 
    x = "",
    title = "Individual Body Mass Distribution (1g - Max)")



# Plot next to bottom temperatures
wigley_b_dist_plot

This makes me lean towards belief that the size distribution is more responsive on shorter time scales to the ambient water temperatures in order to align so well with seasonal temperature swings. To me this suggests that as a whole the community is likely better equipped to respond to the long-term trends if they track larger absolute swings within the year.

Code
wigley_bmspectra_model_df %>% 
  ggplot(aes(bot_temp, b)) +
  geom_point(aes(color = season, shape = survey_area)) +
  geom_smooth(method = "gam", color = "black") +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(title = "Broad Relationship to Ambient Bottom Temperature")

Total Landings

The landings information is the newer GARFO data obtained from Andy Beet. It was delivered as containing only finfish, I have removed freshwater species and blue crabs. It should largely reflect marine finfish species at this point.

Code
wigley_bmspectra_model_df %>% 
  distinct(est_year, survey_area, total_weight_lb) %>% 
  ggplot(aes(est_year, total_weight_lb/1e6)) +
  geom_col(alpha = 0.5) +
  geom_ma(aes(color = "5-Year Smooth"), size = 1, n = 5, ma_fun = SMA, linetype = 1) +
  scale_color_manual(values = "orange") +
  facet_wrap(~survey_area, ncol = 1) +
  scale_y_continuous(labels = label_number(suffix = "M")) +
  labs(y = "Total Landings (live lb.)", x = "Year", color = "Moving Average")

Zooplankton Community Indices

The ecodata package has changed the zooplankton indices around that it carries. It now contains “regime” metrics for 26 of the main zooplankton taxa.

Code
zoo_regimes <- ecodata::zoo_regime
# zoo_anoms <- ecodata::zoo_abundance_anom
zoo_regimes %>% 
  filter(str_detect(Var, "calfin|chaeto|pseudo")) %>% 
ggplot(aes(Time, Value, color = Var)) +
  geom_point(
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(linetype = "5-Year Smooth"),
    size = 1, n = 5, ma_fun = SMA) +
  facet_grid(EPU~.) +
  labs()


Driver Correlations

Driver Modeling

This section is likely to change

The following model has been used to determine the relationship between landings and seasonal bottom temperature on b.

lmer(
  b ~ survey_area*season*scale(roll_temp) + survey_area * scale(log(total_weight_lb)) + (1|est_year), 
  data = wtb_model_df)
Code
#### Model 2: Temperature & Fishing ####


# Drop NA's
wtb_model_df <- drop_na(wigley_bmspectra_model_df, total_weight_lb, bot_temp, b) %>%
  rename(area = survey_area) %>% 
  left_join(area_df) %>% 
  # Do rolling BT within a season * region
  group_by(survey_area, season) %>%
  mutate(
    roll_temp = zoo::rollapply(bot_temp, 5, mean, na.rm = T, align = "right",  fill = NA)) %>% 
  ungroup() %>% 
  mutate(
    yr_num = as.numeric(est_year),
    yr_fac = factor(est_year),
    survey_area = factor(survey_area, levels = area_levels),
    season = factor(season, levels = c("Spring", "Fall")),
    landings = total_weight_lb,
    yr_seas = str_c(season, est_year))


# Make the model
mod2 <- lmer(
  b ~ survey_area*season*scale(roll_temp) + survey_area*scale(log(total_weight_lb)) + (1|est_year), 
  data = wtb_model_df)


# summary(mod2)
# check_model(mod2)

Temperature Relationship

Code
# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    mod2, 
    terms = ~ roll_temp*survey_area*season))   %>% 
  mutate(
    survey_area = factor(group, levels = area_levels),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####
trend_df <- emtrends(
  object = mod2,
  ~survey_area * season,
  var = "roll_temp")


# Drop effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = mod2,
  specs =  ~ survey_area*season,
  var = "roll_temp") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))




# Limit temperature plotting range
actual_range <- wtb_model_df %>% 
  group_by(season, survey_area) %>% 
  summarise(min_temp = min(bot_temp)-2,
            max_temp = max(bot_temp)+2,
            .groups = "drop")



# Plot over observed data
mod2_preds %>% 
  left_join(actual_range) %>% 
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>% 
  left_join(select(mod2_emtrend, survey_area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = survey_area), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wtb_model_df,
    aes(roll_temp, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(survey_area~., scales = "free") +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Exponent of ISD",
       title = "Wigley Species, Body Mass Spectra",
       x = "5-Year Rolling Average Bottom Temperature")

Landings Relationship

Code
# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    mod2, 
    terms = ~ total_weight_lb * survey_area * season))   %>% 
  mutate(
    survey_area = factor(group, levels = area_levels),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####
trend_df <- emtrends(
  object = mod2,
  ~survey_area * season,
  var = "total_weight_lb")
#summary(trend_df, infer= c(T,T))


# Drop effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = mod2,
  specs =  ~ survey_area*season,
  var = "log(total_weight_lb)") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T),
    group = str_c(survey_area, "X", season))




# Plot over observed data
mod2_preds %>% 
  left_join(select(mod2_emtrend, survey_area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wtb_model_df,
    aes(total_weight_lb, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(survey_area~., scales = "free") +
  scale_color_gmri() +
  scale_x_continuous(transform = "log10", labels = label_log(base = 10)) +
  labs(y = "Body Mass Spectra Slope (b)",
       title = "Wigley Species, Body Mass Spectra",
       x = "Total Landings (lb.)")

Driver Model Summary

Code
# mod2 %>% 
#   gtsummary::tbl_regression()   %>% 
#   add_global_p(keep = T) %>% 
#   bold_p(t = 0.05)  %>%
#   bold_labels() %>%
#   italicize_levels() %>% 
#   as_gt() %>% 
#   tab_header(title = "Wigley Community - Bodymass Distribution & Drivers")

Side Stories

Everything below is digging into different nuances of the approach or the community composition. These things are of secondary importance and will be moved out.

Spiny Dogfish Sensitivity

Spiny dogfish constitute a large majority fraction of biomass in the MAB during the Spring survey. Their numbers have also been increasing over time. This acts to shallow the spectra as they are one of the larger species routinely sampled by the trawl survey.

If we remove them from the estimation and re-run the slopes this is what changes.

How much can one species shift the spectra, and is this a definitive/conclusive proof that its the dogfish shallowing the slope in MAB?

Code
# # Run MAB with and without spiny dogfish
# # Run the year, season, area fits for the filtered data
# no_dogs_bmspectra <- group_binspecies_spectra(
#   ss_input = filter(wigley_in, comname != "spiny dogfish"),
#   grouping_vars = c("est_year", "season", "survey_area"),
#   abundance_vals = "numlen_adj",
#   length_vals = "length_cm",
#   use_weight = TRUE,
#   isd_xmin = 1,
#   global_min = TRUE,
#   isd_xmax = NULL,
#   global_max = FALSE,
#   bin_width = 1,
#   vdiff = 2)

# # Save those
# write_csv(
#   no_dogs_bmspectra,
#   here::here("Data/processed/wigley_species_nodogfish_bmspectra.csv"))
Code
# Read them in, do some plotting
no_dogs_bmspectra <- read_csv(
  here::here("Data/processed/wigley_species_nodogfish_bmspectra.csv")) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels))


# Format a little
no_dogs_bmspectra <- no_dogs_bmspectra %>% 
  mutate(est_year = as.numeric(est_year))


# Join them together
dfish_results <- bind_rows(
   list(
     "Wigley Full" = wigley_bmspectra,
     "Wigley w/o Spiny Dogfish" =  no_dogs_bmspectra), 
   .id = "community") %>% 
  mutate(metric = "bodymass_spectra") %>% 
  mutate(
    survey_area = fct_relevel(survey_area, area_levels))




# Get trends for no dogfish
nodfish_mod <- lmer(
  formula = b ~ est_year * survey_area * season + (1|est_year),
  data = no_dogs_bmspectra)
bmspectra_mod_wig <- lmer(
  formula = b ~ est_year * survey_area * season + (1|est_year),
  data = wigley_bmspectra)

# Put predictions
dogfish_sens_predictions <- bind_rows(list(
  "bodymass_spectra-Wigley Full" = get_preds_and_trendsignif(bmspectra_mod_wig),
  "bodymass_spectra-Wigley w/o Spiny Dogfish" = get_preds_and_trendsignif(nodfish_mod)), 
  .id = "model_id") %>% 
  separate(model_id, into = c("metric", "community"), sep = "-") %>% 
  mutate(metric = fct_rev(metric))



# Plot the differences
ggplot() +
  geom_ribbon(
    data = filter(dogfish_sens_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    data = dfish_results,
    aes(est_year, b, color = season),
    size = 0.8, alpha = 0.5) +
  geom_line(
    data = filter(dogfish_sens_predictions, non_zero == TRUE),
    aes(x, predicted, color = season),
    linewidth = 1, alpha = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(survey_area ~ metric*community, scales = "free") +
  labs(
    x = "Year",
    y = "Exponent of Size Spectra",
    title = "Shift in Spectra Trend(s) with Exclusion of Spiny Dogfish")

Story Thoughts

Below is a table of some papers operating in the same area or around the same topic that are particularly relevant:

Code
shelf_papers <- tribble(
  ~"Author", ~"Year", ~"Conjecture",
  "Friedland et al.", 2024, "Body size has declined despite more abundance & biomass. Competition is preventing larger sizes.",
  "Krumsick", 2020, "Empirical slopes from trawl survey in labrador sea were steeper than theoretical from nitrogen stable isotope ratios"
)


shelf_papers %>% gt()
Author Year Conjecture
Friedland et al. 2024 Body size has declined despite more abundance & biomass. Competition is preventing larger sizes.
Krumsick 2020 Empirical slopes from trawl survey in labrador sea were steeper than theoretical from nitrogen stable isotope ratios

There are a couple core ideas about what changes have occurred or have been ocurring over the Northeast Shelf with relevance to the community that is sampled by the trawl survey: